library(readxl)
library(dplyr)
library(stringr)
library(ggplot2)

WEINKLEDATA <- "data/41893_2018_165_MOESM2_ESM.xlsx"

dir.create("out", showWarnings = FALSE)
dir.create("out/figures", showWarnings = FALSE)

theme_set(theme_minimal(base_size = 10))

# Data Cleaning:

ALLSTORMS <- read_excel(WEINKLEDATA, sheet = "Storm Damages")[, 1:8]
names(ALLSTORMS) <- c("year", "storm", "name", "cat", "states", "basedamage", "pl18", "cl18")

ALLSTORMS <- ALLSTORMS %>%
  filter(!is.na(year)) %>%
  mutate(year = as.integer(year),
    category = as.integer(cat),
    basedamage = as.numeric(basedamage),
    pl18 = as.numeric(pl18),
    cl18 = as.numeric(cl18),
    states = str_squish(states),
    name = str_squish(name),
    imputed = is.na(basedamage))

imputedconstants <- ALLSTORMS %>%
  filter(imputed) %>%
  count(category, pl18) %>%
  arrange(category, pl18)
print(imputedconstants)

sev <- ALLSTORMS %>% filter(!imputed)

# Some Checks

stopifnot(nrow(ALLSTORMS) == 197,
  nrow(sev) == 160,
  abs(min(ALLSTORMS$pl18) - 2785748) < 1)

# Fitting the distribution for severity

x  <- sev$pl18 / 1e9

# Note that this changes units to be in Billions ($USD).

lx <- log(x)
mu <- mean(lx)
sigma <- sd(lx)
EX  <- exp(mu+sigma^2/2)
EX2 <- exp(2*mu+2*sigma^2)
SDX <- sqrt((exp(sigma^2)-1)*exp(2*mu+sigma^2))

cat("n:", length(x), ".\n")
cat("mu:", round(mu, 4), ".\n")
cat("sigma:", round(sigma, 4), "\n")
cat("E[X]: $", round(EX, 3), "billion.\n")
cat("SD[X]: $", round(SDX, 3), "billion.\n")
cat("median: $", round(exp(mu), 3), "billion.\n")
cat("sample mean: $", round(mean(x), 3), "billion.\n")


# Testing whether removing imputed storms makes a difference

lxall <- log(ALLSTORMS$pl18 / 1e9)
effectofimputed <- data.frame(sample = c("Only Include Observations", "All Included"),
  n = c(length(lx), length(lxall)),
  mu = round(c(mean(lx), mean(lxall)), 4),
  sigma = round(c(sd(lx),   sd(lxall)),   4))
print(effectofimputed, row.names = FALSE)

# Testing the fit of the lognormal distribution

n <- length(x)

qq <- data.frame(empirical=sort(x), theoretical=qlnorm((seq_len(n)-0.5)/n, mu, sigma))

pqq <- ggplot(qq, aes(theoretical, empirical))+geom_abline(slope = 1, colour = "firebrick")+
  geom_point(size = 1.4, alpha = 0.7)+scale_x_log10() + scale_y_log10()+
  labs(x = "Quantiles", y = "Observed Loss ($ Billions)", title = "Q-Q Plot for the Fit of Lognormal to Losses from US Hurricanes")

ggsave("out/figures/fig_qq_lognormal.pdf", pqq, width = 6, height = 4.5)

# Survival Plot

survival <- data.frame(
  loss = sort(x, decreasing = TRUE),
  survival = seq_len(n) / (n + 1)
)

psurvival <- ggplot(survival, aes(loss, survival))+ geom_point(size = 1.3, alpha = 0.7)+
  stat_function(fun = function(q) plnorm(q, mu, sigma, lower.tail = FALSE),
  colour = "firebrick")+scale_x_log10() + scale_y_log10()+
  labs(x = "Normalized Loss ($ Billion)", y = "P(X > x)",
  title = "An Empirical survival Function versus the Fitted Lognormal")

ggsave("out/figures/fig_survival.pdf", psurvival, width = 6, height = 4.5)


# Simulation of Null Envelope 

set.seed(67420)
simulated <- replicate(10000, sort(rlnorm(n, mu, sigma), decreasing = TRUE))
observed <- sort(x, decreasing = TRUE)

band <- data.frame(survival=seq_len(n)/(n + 1),
  observed = observed,
  low = apply(simulated, 1, quantile, 0.025),
  high = apply(simulated, 1, quantile, 0.975))
band$outside <- band$observed < band$low | band$observed > band$high

pband <- ggplot(band, aes(y = survival)) +
  geom_ribbon(aes(xmin = low, xmax = high), fill = "grey85")+geom_point(aes(x = observed), size = 1.3)+
  scale_x_log10()+scale_y_log10()+
  labs(x = "Normalized Loss ($ Billion)", y = "P(X > x)",
  title = "Observed Losses versus a Lognormal Null Envelope",
  subtitle = "The Gray Band Represents a 95% Range from 10000 Simulated Samples of n=160")

ggsave("out/figures/fig_null_envelope.pdf", pband, width=6, height =4.5)

cat("Of those Outside the 95% band, the Top 10:", sum(band$outside[1:10]),
    "Top 30:", sum(band$outside[1:30]),
    "All:", sum(band$outside), "of", n, "\n")

cat("The Band for the Largest Loss: [$", round(band$low[1], 1),
    "Billions, $", round(band$high[1], 1), "bn], observed $", round(band$observed[1], 1), "bn\n")


# Testing if/why normalization matters:

methodparams <- data.frame(method = c("PL18", "CL18"),
  mu = round(c(mean(log(sev$pl18/1e9)), mean(log(sev$cl18/1e9))), 4),
  sigma = round(c(sd(log(sev$pl18/1e9)),   sd(log(sev$cl18/1e9))),   4)
)
print(methodparams, row.names = FALSE)


# The combined model

lambda <- 1.663
knb <- 7.358


pmf <- data.frame(
  n = 0:6,
  poisson = round(dpois(0:6, lambda), 4),
  negbin = round(dnbinom(0:6, size = knb, mu = lambda), 4)
)
print(pmf, row.names = FALSE)
cat("P(N>=5) Poisson:", round(1 - ppois(4, lambda), 4),
    "Negative Binomial:", round(1 - pnbinom(4, size = knb, mu = lambda), 4), "\n")

ES <- lambda*EX
varSpois <- lambda*EX2
varSnb <- lambda*EX2+(lambda^2 / knb)*EX^2

cat("E[S]: $", round(ES, 3), "bn\n")
cat("Var[S] Poisson:", round(varSpois, 1), "\n")
cat("Var[S] NB:", round(varSnb, 1), "\n")
cat("variances Ratio:", round(varSnb / varSpois, 5), "\n")

set.seed(67420)
NSIM <- 1000000

sim <- function(N) {vapply(N, function(n) if (n == 0) 0 else sum(rlnorm(n, mu, sigma)), numeric(1))}

Spois <- sim(rpois(NSIM,   lambda))
Snb <- sim(rnbinom(NSIM, size = knb, mu = lambda))

cat("Simulated E[S]:", round(mean(Spois), 3),
    "E[S] Given by Formula:", round(ES, 3), "\n")
cat("Simulated Var[S]:", round(var(Spois), 1),
    "Var[S] Given by Formula:", round(varSpois, 1), "\n")
cat("Simulated P(S=0):", round(mean(Spois == 0), 4),
    "exp(-lambda):", round(exp(-lambda), 4), "\n")


summariseS <- function(S, label) {
  data.frame(model = label,
    AAL = mean(S),
    AEP100 = as.numeric(quantile(S, 0.99)),
    AEP250 = as.numeric(quantile(S, 0.996))
  )
}

results <- rbind(summariseS(Spois, "Poisson"), summariseS(Snb,   "Negative binomial"))
results$delta_AAL_pct <- 100 * (results$AAL / results$AAL[1] - 1)
results$delta_250_pct <- 100 * (results$AEP250 / results$AEP250[1] - 1)

print(cbind(model = results$model, round(results[-1], 3)), row.names = FALSE)

set.seed(67420)
NBOOT <- 1000

aep250draws <- vapply(seq_len(NBOOT), function(j) {
  lb <- log(sample(x, replace = TRUE))
  m  <- mean(lb); s <- sd(lb)
  N  <- rpois(2000, lambda)
  S  <- vapply(N, function(k) if (k == 0) 0 else sum(rlnorm(k, m, s)), numeric(1))
  as.numeric(quantile(S, 0.996))}, numeric(1))

cat("Point Estimate: $", round(results$AEP250[1], 2), "bn\n")
cat("95% Interval  : [$", round(quantile(aep250draws, 0.025), 2),
    "Billion, $", round(quantile(aep250draws, 0.975), 2), "bn]\n")
cat("Width as a % of the Point Estimate:",
    round(100 * (quantile(aep250draws, 0.975) - quantile(aep250draws, 0.025))/results$AEP250[1], 1), "%\n")
cat("compare to the Poisson-vs-NB difference of",
    round(results$delta_250_pct[2], 2), "%\n")

set.seed(67420)
bootpois <- replicate(500, quantile(sample(Spois, replace = TRUE), 0.996))
bootnb <- replicate(500, quantile(sample(Snb, replace = TRUE), 0.996))
print(quantile(bootpois, c(0.025, 0.975)))
print(quantile(bootnb, c(0.025, 0.975)))